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A density functional theory for a macroion suspension is examined, where the excess free energy 
corresponds to the macroion self energy arising from the polarisation of the supporting electrolyte 
solution. This is treated within a linearised or Debye-Hiickel approximation. The model predicts 
liquid-liquid phase separation at low ionic strength. The interface structure and surface tension 
between coexisting phases is calculated using a variational approximation. Results are also obtained 
for structure factors, which are shown to obey the Stillinger-Lovett moment conditions. As one 
approaches the critical points, the structure factors may diverge at a non-zero wavevector, indicating 
that the critical points could be replaced by charge-density-wave phases. 



I. INTRODUCTION 

The phase behaviour of charged colloidal suspensions 
at low ionic strength has attracted much experimental 
and theoretical interest. Observations of void structures 
and other phenomena Q, 0] motivated a number of theo- 
retical studies which attributed the anomalous behaviour 
to phase separation between colloid-rich and colloid-poor 
phases HHH®. 

Several reviews are available M, |8l . 
The original theoretical explanations have come under 
strong attack for using a Debye-Hiickel linearisation ap- 
proximation which is, at best, at the margin of its va- 
lidity. Various attempts to patch this up have left the 
situation unclear. Cell model calculations using Poisson- 
Boltzmann theory indicate the original predictions are an 
artefact of the linearisation approximation 0, 0] . The 
Debye-Hiickel approximation can be improved by tak- 
ing into account counterion condensation, in which case 
the phase transition may or may not be recovered de- 
pending on the approximation scheme used Other 
approaches such as extended Debye-Hiickel theory \\% . 
symmetrised Poisson-Boltzmann theory 0, 0], 'boot- 
strap' Poisson-Boltzmann theory [l^, and a systematic 
expansion into two- and three-body interactions |l6l Il7| , 
all indicate that a phase transition can occur, as do sev- 
eral integral equation studies |l8l [19| . The experimental 
situation is also uncertain since a plausible alternative 
explanation has been suggested y|, in which the voids 
correspond to regions occupied by dilute, highly extended 
(and therefore effectively invisble) polyelectrolyte chains 
which have been shed by the latex colloids. 

Simulation methods struggle to approach these prob- 
lems because of the disparity in size between the 
macroions and the small ions, and the need to handle 
the electrostatic interactions. Nevertheless, convincing 
evidence has been found for liquid-liquid phase separa- 
tion in a macroion system at lower dimensionless tem- 
peratures [23, COI ■ Experimentally this corresponds to a 
solvent with a lower dielectric constant than water (but 
one in which the ions still disperse). Charged colloids in 
such solvents exhibit many interesting phenomena • 

Thus, whilst the weight of evidence perhaps suggests 
that aqueous charge-stabilised colloidal suspensions may 



not show genuine liquid-liquid phase coexistence, it is 
absolutely clear that there will be phase coexistence be- 
tween condensed and dilute colloidal phases at small 
enough dimensionless temperatures. In this sense, the 
problem resembles the much-studied restricted primitive 
model (RPMj, whose phase behaviour is now well estab- 
lished |M El III S3- 

In situations where genuine phase coexistence obtains, 
one can go on to ask questions about the surface tension 
and electrical structure of the interface between the co- 
existing phases. Answers to these questions may prompt 
new avenues for experimental investigation of real sys- 
tems. Previously, Knott and Ford compute the surface 
tension using square-gradient theory, but discard the pos- 
sible electrical structure at the interface [2]j . The present 
work approaches this problem within the context of a 
density functional theory, motivated by the earlier study 
in Ref. p| (see also Appendix A). It places the phe- 
nomenological remarks made in this earlier work on a 
sounder footing. The analysis in Ref. [j| suggests that 
the macroion self energy is the dominant contribution to 
the excess free energy, similar to an early insight by Lang- 
muir (2^. In the present work therefore, the rather gross 
simplification has been adopted in which the macroion 
self energy is the only contribution to the excess free en- 
ergy. Moreover this self energy is computed in a simple 
closed form using Debye-Hiickel theory, and is thus also 
based on the much-criticised linearisation approximation. 
Nevertheless I argue that it is instructive to proceed, be- 
cause of the rich phenomenology that is revealed. 

The model predicts phase separation at low dimension- 
less temperatures and low ionic strengths, and in quan- 
titative terms stands reasonable comparison with some 
of the other approaches. The physics of the phase sep- 
aration lies in the dependence of the macroion self en- 
ergy on the local ionic strength: macroions drift towards 
regions of high ionic strength, which by charge neutral- 
ity are regions where other macroions have also congre- 
gated. Within the linearisation approximation, the ef- 
fect grows without bound as the macroion charge is in- 
creased, and thus the mechanism can drive phase sepa- 
ration at sufficiently large macroion charges. In reality, 
non-linear effects (counterion condensation) limit the ef- 
fective macroion charge |29|. and therefore this mecha- 
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nism is probably insufficient in itself to drive phase sep- 
aration in real systems. Undoubtably though it is still a 
contributing factor, operating in conjunction with other 
effects such as correlated fluctuations in the counterion 
clouds around macroions and the sharing of counterions 
between macroions 0, H0| ■ 

The model is constructed in the form of a density func- 
tional theory. Thus, as well as making predictions for 
phase separation, it can be used to solve for the den- 
sity profiles and the surface tension between coexisting 
phases. The results obtained here are in accord with 
typical expectations for soft condensed matter systems 
|3l| . and were summarised in an earlier publication |32j . 
In addition, I also discuss the predictions that the the- 
ory makes for the structure factors. These are found to 
obey the Stillinger-Lovett moment conditions [33l l34| . 
although it turns out this is not a stringent test of the 
theory. Intriguingly, I find that the structure factors may 
diverge at a non-zero wavevector as one approaches the 
critical points. This suggests the possibility that the crit- 
ical points in these s yste ms may be replaced by charge- 
density- wave phases |35| . This phenomenological possi- 
bility in charged systems was first suggested by Nabu- 
tovskii and coworkers I23l 13a, 1371 . 



II. SPECIFICATION OF THE MODEL 

The underlying model of the macroion suspension used 
here is a primitive model commonly deployed for this 
kind of problem. The 'primitive' aspect is that the sol- 
vent is treated as a structureless dielectric continuum in 
which the macroions and small salt ions are embedded. 
The macroions are treated as spheres of (positive) charge 
Z, diameter a, and number density p rn (volume fraction 

= ira 3 p m /6). The salt ions are univalent counterions 
and coions at number densities p- and p+ respectively. 

1 suppose there is only one species of counterion. The 
size of the salt ions is assumed to be small enough to be 
irrelevant. The dielectric continuum is characterised by 
a Bjerrum length Ib so that the electrostatic interaction 
energy between a pair of univalent charges separated by 
a distance r is Ib/i", in units of k^T where fee is Boltz- 
mann's constant and T is the temperature. For water 
at room temperature, Ib ~ 0.72 nm. The model is com- 
pletely parametrised by the dimensionless ratio a /Is and 
the charge Z. It is often convenient to pretend that the 
dielectric permittivity of the background is independent 
of temperature, in which case Ib ~ 1/T. This means that 
a /Ib can be regarded as a dimensionless temperature. 

The density functional theory is specified by giving the 
free energy F as a functional of the spatially varying 
number densities p m ( r ) an d P±( r ) US- The functional is 
decomposed into ideal, mean-field, and correlation con- 
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The first term is the ideal term: e is the base of nat- 
ural logarithms and the pf are unimportant base units 
of concentration related to the definition of the standard 
state [33. The second term is a mean-field electrostat- 
ics term: p z (r) = J2i z iPi( r ) * s the local charge density 
with Zi = {Z, 1, —1} as i = {771, +, — }, and a factor 1/2 
allows for double counting. The third term (correlation 
term) represents the excess free energy. As discussed 
above, only the macroion self energy f m is included in 
this term. This is computed using Debye-Huckel theory 

11111111113, 
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where re(r) is a local inverse Debye screening length. This 
is defined in terms of an average local ionic strength, 
Pj(r), through 



[«(r)] 2 = 8^ B P/(r), 
p I (r)=fd 3 r'w(\r-r'\)p I (r') 
Pl (r') = [p + (r') + p_(r')}/2. 



(3) 



The ionic strength includes the counterions and salt ions, 
but not the macroions. In principle, allowance should be 
made for the macroion excluded volume, but this effect 
is of secondary importance and for simplicity has been 
omitted. 

The smoothing kernel in the second of Eqs. J3J is nor- 
malised so that / d 3 rw(r) = 1. Here I use 



w(r) = (nacr 2 ) 3 / 2 exp[-r 2 /(cw 2 )]. 



(4) 



This is an arbitrarily chosen function 43], of range 
The argument below suggests that the parameter 
a should be of order unity and for the most part I will set 
a = 1 in the calculations. Eqs. Q-0 completely specify 
the density functional theory, and everything discussed 
below can be derived from them. 

The decomposition into ideal, mean field, and correla- 
tion contributions is a standard approach pfl l45l l4t| . 
The approximation made for the correlation term de- 
serves more discussion though. The only piece of physics 
that has been incorporated is the macroion self energy. 
This has a non-trivial dependence on the local ionic 
strength since each macroion polarises the surrounding 
electrolyte and becomes surrounded by a 'double layer'. 



3 



This dependence causes macroions to drift towards re- 
gions of high ionic strength, as discussed already. 

The physical reason for introducing a smoothing kernel 
is that one can derive the self energy by integrating out 
the small ion degrees of freedom, with the main contribu- 
tion coming from variations on length scales correspond- 
ing to the structure in the double layer Thus only 
variations in ionic strength on length scales > a should 
be included in the model. The smoothing kernel is a de- 
vice for achieving this. This argument also motivates the 
choice for a in Eq. 

In section Ivl below, it is found that the theory is not 
well behaved if one uses a 'point model' where the depen- 
dence is on the ionic strength at, say, the centre of the 
macroion (the first of Eqs. Q with ~p~j replaced by pi). 
This provides a second technical reason to make the self 
energy depend on a smeared ionic strength. 

The potential energy of a small ion at the surface of 
the macroion, in units of fceT, is ±ZIb/ct. Eq. J5J uses 
the Debye-Hiickel expression for the self energy, which 
assumes ZI-q/ct <C 1. The expression becomes increas- 
ingly inaccurate for Zl^/a > 1, and its use has been 
the subject of strong criticism as discussed above. Since 
the interesting effects are found only at larger values of 
ZIb/ct, one should interpret the quantitative results with 
caution. 



III. BULK PHASE BEHAVIOUR 
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In this section, I shall consider the bulk phase be- 
haviour predicted by the free energy of Eqs. Q-Q- This 
is a homogeneous situation in which the density variables 
lose their spatial dependence. In this limit, one can prove 
that the mean field term should be replaced by a condi- 
tion of bulk charge neutrality, p z = J2i ZiPi = [T3,E3l- 

The required charge neutrality condition can be im- 
posed in two ways. The first route is to add a term 
^k^T J2 i ZiPi to the free energy, where ipk^T is a La- 
grange multiplier. This approach has the advantage of 
making a close connection to the density functional the- 
ory. Taking this approach, the free energy becomes 
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where V is the system volume and k 2 — 4itIb{p+ + /?-)■ 
The distinction between the smoothed and unsmoothed 
ionic strength disappears in the homogeneous limit. In 
this approach the pi are treated as three independent 
density variables. At the end of any calculations, i/j is 
adjusted to get J2i z %Pi — 0- The value of tp depends on 
the state point under consideration. 

The second way to enforce charge neutrality is to elim- 
inate one of the density variables. Since this is numer- 
ically quite convenient, it is the approach that shall be 
adopted in the rest of this section. At this point one 
can recognise that the coions come from added salt and 



FIG. 1: (a) Universal phase behaviour in the absence of salt, 
predicted by Eq. @ (upper curve) . The simulation results of 
Rescic and Linse 1211 are also shown (lower curve with marked 
points) . (b) Behaviour of the critical point at Z = 10 3 as salt 
is added. The dashed line corresponds to the parameters used 
in Fig. |3 



write p- = Zp m + p s and p+ = p s , where p s is the added 
salt concentration. The free energy is given by Eq. (J5J 
but with -0 = 0, and p± substituted by the above ex- 
pressions. There are now only two independent density 
variables and the phase behaviour can be represented in 
the (p m ,p s ) plane. 

I now discuss the phase behaviour predicted by this 
free energy. Firstly, in the absence of salt some additional 
simplifications can be made. In the limit p s — > 0, the free 
energy can be written in a dimensionless form as 
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where (j) IS the macroion volume fraction. To get to this 
point, I have assumed that 2>1 and hidden some con- 
stants and terms strictly proportional to p m since they 
do not affect the phase behaviour. 

Eq. JfJJl predicts the dependence on a/l-B and Z is 
through the single combination Zl^/a (there is no rea- 
son to suppose that this should be the case in a more 
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FIG. 2: Phase behaviour at Z = 10 3 , a = 100 nm and 
Ib = 0.72 nm, corresponding to the dashed line in Fig. Q 
The miscibility gap is bounded above and below by critical 
points. The dashed tie line is the one for which the interfacial 
properties are reported in Figs. [2HS1 



accurate theory). This is the same parameter that quan- 
tifies the accuracy of the Debye-Hiickel linearisation ap- 
proximation. The inverse of this, a/^Zl-g), is propor- 
tional to the dimensionless temperature discussed above. 
Fig. ^a) shows the universal phase behaviour predicted 
by Eq. © as a function of the macroion volume frac- 
tion and ajiZl-Q). At small enough values of <t/(ZIb), 
a two phase region is encountered in the phase diagram. 
The two phase region corresponds to phase coexistence 
between macroion rich and macroion poor phases. The 
identities of these phases merge at a critical point located 
at (j) « 9.18 x 10~ 3 and a/(Zl B ) « 0.132. 

One can compare this with the simulation results of 
Rescic and Linse for Z = 10 macroions [2l|. They also 
find a two phase region on lowering temperature, with a 
critical point located at 4> w 0.17 and <j/(ZIb) ~ 0.077. 
Whilst the phenomenology is the same, the numerical val- 
ues are somewhat different from the prediction of Eq. 10. 
Not unexpectedly, the present model is too crude to ob- 
tain quantitatively reliable results. An analogy can be 
made with the application of Debye-Hiickel theory to the 
restricted primitive model (RPM) [23l I2M lisf . In this 
case too, Debye-Hiickel theory correctly suggests a re- 
gion of phase separation at low temperatures but errs 
in terms of quantitative predictions. Interestingly, in 
terms of accuracy of prediction, the present theory is not 
much worse than symmetrised Poisson-Boltzmann theory 
or the mean spherical approximation |13l Il9l |49| . 

I now turn the effect of added salt, and analyse the 
predictions of the full free energy in Eq. JSJ. In general, 
as salt is added, the critical point in Fig.^a) first moves 
to higher dimensionless temperatures, passes through a 
maximum, and then starts to move to lower dimension- 
less temperatures again. This non-monotonic behaviour 
is shown in Fig. db) for Z = 10 3 . A similar effect 
of added salt is seen in a number of other approaches 




FIG. 3: Macroion volume fraction (top) and small ion con- 
centrations (bottom) through the interface corresponding to 
the dashed tie line in Fig. [5] 



I40L I50L l&lj, |&2j. In the presence of added salt, it is 
no longer true that the dependence on Z and a /Ib can 
be combined into a single parameter, however for com- 
parison with the phase behaviour in the absence of salt, 
Fig. ^b) shows the behavior as a function of aj{ZlB) at 
this fixed value of Z . 

The re-entrant behaviour means that for parame- 
ters such as those corresponding to the dashed line in 
Fig- ffib), there are two critical points in the (p m ,p s ) 
plane, and one encounters a re-entrant single phase re- 
gion at low added salt. The dashed line in Fig.^b) is for 
Z = 10 3 , (j = 100 nm and l B = 0.72 nm, and the corre- 
sponding phase behaviour in the {p m , Ps) plane is shown 
in Fig. |21 It is seen that the two phase region appears as 
a miscibility gap in this representation. 

As <j/Ib is increased or Z is decreased, the two critical 
points move towards each other and finally disappear at 
a double critical point, or hypercritical point |53| . For ex- 
ample, for Z — 10 3 the double critical point corresponds 
to the maximum of the solid line in Fig. ^b), where 
<j/(Zl B ) « 0.145, <j> w 1.04 x 10~ 2 and p s w 8.98 pM. 

The bulk phase behaviour predicted by Eq. thus 
closely resembles that predicted by various other ap- 
proaches, including the theory discussed in Ref. |5|. 
Many approaches, including the present one, do not con- 
sider the formation of ordered phases (colloidal crystals). 
These can arise from the strong macroion-macroion in- 
teractions. The possibility of ordered phases has been 
considered by van Roij and coworkers 0, [13, E3 though. 
They find that ordered phases can appear in the vicin- 
ity of the miscibility gap in which case a richer phase 
behaviour can result. 
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FIG. 5: Excess grand potential density (solid line) and elec- 
trostatic component thereof (dashed line) corresponding to 
the ion density profiles shown in Fig. 03 



FIG. 4: Charge density (top), electric field (middle), and elec- 
trostatic potential (bottom) corresponding to the ion density 
profiles shown in Fig. 03 



IV. INTERFACIAL PROPERTIES 

A major use of the density functional theory in the 
present context is to calculate the macroion and small 
ion density profiles through the interface between two 
coexisting phases, and to compute the surface tension. In 
order to set the problem up, it is convenient to introduce 
the grand potential |3S| 



n = F- Jd 3 r Mi( 



(7) 



where pi are the chemical potentials of the three species, 
and F is defined in Eqs. At this point it is also 

convenient to rewrite the mean field term in Eq. l(T|l. De- 
fine a dimensionless electrostatic potential 



1>(r) = l B J d 3 r 
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so that the mean field term in Eq. can be written 

l fJd^r^-f^= 1 -J d ^ { r )Pz{ r). (9) 

By direct substitution, one verifies that the potential de- 
fined by Eq. JSJ solves the Poisson equation 



V 2 ip + 4:irl B p z = 0. 



(10) 



Using this and Green's first identity [5J|, the mean field 
term now becomes 

l -Jd?v^{v) Pz {v) = ^ I -Jd?v\V^\ 2 . (11) 



This is recognised as the electric field energy since Vtjj is 
essentially the electric field strength. One can now define 
a grand potential density w(r) such that VL = fd 3 ru(r) 
and 
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where the explicit dependence on the spatial co-ordinate 
has been suppressed. For a homogeneous system, to = —p 
where p is the pressure. 

Setting Sfl/Spi(v) = and using Eq. (JSJ) gives 
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In principle, these non-linear integral equations can be 
solved to find the ion density profiles. Here a variational 
approximation has been adopted in which f2 is minimised 
with respect to parameters in trial functions which spec- 
ify the ion density profiles. More details of the numerical 
approach are given in Appendix B. 

I now suppose that all the variation occurs in one di- 
rection x normal to the interface. At large distances 
from the interface, x — > ±oo, the number densities ap- 
proach those corresponding to the coexisting bulk phases. 
The grand potential density approaches a constant value 
w(±oo) equal to (minus) the pressure, and therefore the 
same in coexisting phases. The surface tension 7 can 
therefore be identified as the excess grand potential per 
unit area 

y = f^dxlu -u>(±oo)]. (14) 
The chemical potentials derived from Eq. (JSJ are 

Mi 
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Comparison with Eq. (|13fl shows that ip in this expres- 
sion is simply the limiting value of t/>(r) in the case of 
a homogeneous system [5.5J- For the interface problem, 
one has two limiting values, ip(±oo). The difference 
Aip — ip(oo) — ip{—oo) arises because of the electrical 
structure at the interface. It is a liquid-liquid junction 
potential analogous to the Donnan potential that ap- 
pears across a semi-permeable membrane |56j| . Since tp 
in Eq. Q15[) is determined by the bulk densities, the dif- 
ference Aip can be calculated without having to solve for 
the interface structure. In fact, because of the symmetric 
way that p± enters into the excess free energy, a simple 
expression obtains, 



A^llnf^M.^j- 00 ; 
2 \p_|_(oo) P-\— oo) 



(16) 



This method of calculating the junction potential was 
used in Ref. Q- 

One question remains: what should be used for the 
chemical potentials in these calculations? The sim- 
plest answer is to compute the chemical potentials from 
Eq. JTSJ, setting tp — and using the bulk densities cor- 
responding to either one of the coexisting phases. This 
works because global charge neutrality means Eq. i|14|) 
for the surface tension is unaffected by the value of ip in 
Eq. I|15fl • Hence we are free to set tp = in either of the 
coexisting phases. 

I now turn to the results. Fig. [31 shows representative 
density profiles for the macroion and small ions through 
the interface between the coexisting phases, correspond- 
ing to the highlighted tie line in Fig. [21 The profiles 
interpolate smoothly between the coexisting bulk densi- 
ties. Fig. 01 shows the detailed electrical structure at the 
interface. The upper plot shows that the charge density 
p z = ZpM+p+ — p- has a dipolar structure. Correspond- 
ingly there is a localised electric field, shown in the middle 
plot, and a smooth jump of ~ 20.5 mV in the elec- 
trostatic potential, shown in the lower plot. This is the 
junction potential which can also be calculated directly 
from the coexisting bulk densities as in Eq. Ijlfcifl . This 
electrical structure is in accord with general expecata- 
tions for charged systems 0, • 

Fig. 03 shows the grand potential density and the 
electrostatic component thereof — the second term of 
Eq. 1)1 2|l — as a function of distance through the inter- 
face. For this particular case the area gives 7 « 0.727 x 
{k-oT/a 2 ). The order of magnitude of this should not 
come as a surprise since er and k^T are the only relevant 
length and energy scales in the problem. Inserting ac- 
tual values, 7 ~ 0.3 /iNnr 1 , which is typical for for soft 
matter interfaces |31| 

Fig. [HI shows how the surface tension and interface 
width vary as one approaches the upper critical point 
in Fig. [21 The width d is defined operationally as d 2 = 
(x 2 ) - (x) 2 , where (...)= J^. . .)p(x) dxj JZoP( x ) dx ' 
w\ihp(x) = \oj(x) — w(±oo)| 2 . These results are obtained 
by repeating the calculations underlying Figs. [3HS1 f° r a 
sequence of tie lines approaching the critical point. They 
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FIG. 6: (a) Surface tension as a function of salt chemical po- 
tential, (b) Interface width (solid line) and correlation lengths 
(dashed lines) as a function of salt chemical potential. In 
both, the salt chemical potential is expressed as a normalised 
distance from the upper critical point. 



are reported as a function of the distance from the critical 
point, expressed in terms of a normalised salt chemical 
potential. Fig. Etb) also shows the correlation lengths 
£± in the coexisting phases determined from the expo- 
nential decay of the density profiles into the bulk phases 
(see Appendix B). As the critical point is approached, 
these approach each other, and diverge in the same way 
as the interface width. Fig.[B|reveals that the surface ten- 
sion and length scales are in accord with expected scaling 
behaviour for a mean-field theory |58j|. 

What happens at the lower critical point in Fig. [21 
though? The next section shows that this is a non-trivial 
question with perhaps an unexpected answer. In the cal- 
culations in the current section, I have assumed that the 
interface profiles smoothly interpolate between the coex- 
isting phases. Indeed, this is the basis of the numerical 
method detailed in Appendix B. However, such an ap- 
proach rules out the possibility of oscillatory behaviour 
in the density profiles (or to be precise, the numerical 
methodology is inappropriate for this scenario). At lower 
salt concentrations though, one can enter a region where 
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FIG. 7: (a) Macroion structure factors at (p — 0.03, and p s = 
4pM (A) and 30 fj,M (B). The inset shows the same curves 
in a double- logarithmic plot. The normalisation is such that 
Smm/pm — > 1 as q — > co. (b) Phase diagram augmented 
by the spinodal line (dashed), the Lifshitz line (dotted), and 
the region where the macroion structure factor diverges at a 
non-zero wavevector (shaded). 



oscillatory behaviour is expected. These considerations 
are made mathematically precise in the next section. 



V. STRUCTURE FACTORS 

The structure factors in a homogeneous system can 
be determined from a density functional theory (DFT) 
by functional differentiation [3g. Where accurate struc- 
ture factors are already known, typically from a combina- 
tion of simulation and integral equation approaches, this 
can be used to constrain the DFT. In the present case 
for example, one could try to constrain w(r) in Eq. lj4*)l. 
However accurate structure factors are not known for this 
problem, and furthermore the DFT has been constructed 
to include only the macroion self energy. Thus it does not 
make sense to constrain the DFT and the present section 
simply reports the structure factors that are predicted 
from the theory as given in Eqs. C3>~®- 



= P&j + PiPjhij(q) 



(17) 



where 
fd 3 



r e 



i and 



j run over {m, +,— } and hij(q) = 
(r) is the Fourier transform of the pair cor- 
relation functions hij(r) = gij(r) — 1. Reciprocal space 
quantities will be denoted by a tilde. The bulk densities 
Pi are constants, fixed by the choice of state point. Devi- 
ations away from these will be denoted by Api. Eq. i(T7|) 
uses the normalisation Sij(q) — > piSij as q — > 00, which 



simplifies some of the expressions below 

To obtain the structure factor matrix, start by defining 
the real-space function 



1 



5 2 F 



k^T \8pi(r)Spj(r') / pi(r)^p 



(18) 



where F is the full free energy. The limit of a homoge- 
neous system is taken after the functional differentiation 
step so that S^ 1 only depends on |r — r'| as indicated. 
Transforming to reciprocal space, one can show that 

§r j \q)=Jd 3 re-^ r S^(r) (19) 
is simply the matrix inverse of Sij , 

z2,SiiSjk = S ik- (20) 
j 

These results follow by combining the Ornstein-Zernike 
relation for a multicomponent mixture in reciprocal 



space, 



Ijk 



where cv,- are the di- 



rect correlation functions [59j, with the DFT result that 
Cjj = — (1/ ' kBT)S 2 F cx /5pjSpi where F cx is the excess free 
energy |38l |. 

The route to the structure factors offered by Eqs. I|18l) - 
(|20|l is based on 'classical' arguments jE^- One can also 
make the connection via field theoretical methods. Ex- 
panding the free energy functional to second order gives 



AF 1 



knT 2 



Jd a rd 3 r'J2^Pi(^Pj(r')S^(\r-r'\ 



■-4 21 ) 

where is defined by Eq. (JTHJ- It follows that jSll 

(A Pi (r)A ft (r'))=^-(|r-r'|) (22) 

where <SV, (r) = J rf 3 q/(27r) 3 e lq ' r :(q) is the structure 
factor matrix expressed as a real space quantity. Al- 
though care has to be taken at the point r = r', one can 
easily show that the density-density correlation function 
on the left hand side of Eq. Q22JI is the same as the Fourier 
transform of the right hand side of Eq. Ijl7|l . 

The Stillinger-Lovett moment conditions constrain the 
behaviour of the structure factors in reci pro cal space in a 
particularly clear manner pi III lil Itil . Firstly, the 
zeroth-moment conditions express perfect screening and 
are / d 3 r J2i z iPi9ij( r ) = ~ z j f° r 3 = {«,+>-}■ Using 
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charge neutrality and assuming the structure factors are 
regular at q — 0, one can easily show that this implies 



z i S lJ (q) = 0(q 2 



(23) 



The second-moment condition is 

/ d 3 rr 2 J2ij ZiZjPiPj9tj(r) = -3/(2nl B ). This con- 
strains the long wavelength behaviour of the charge- 
charge structure factor, 



(47ri B ) 



0{q A ) 



(24) 



In real space, this means that (Ap z (r)Ap 2 (r')} ~ 1-q/\y — 
r'| for |r— r'| — * oo. Thus charge density fluctuations van- 
ish with the Coulomb law at large distances, correspond- 
ing to the fact that the electrostatic energy dominates in 
the free energy for long-wavelength density fluctuations 
unless they happen to be charge- neutral |64| . 

I now apply the formalism of Eqs. (jl8|) (|20JI to the 
present DFT defined in Eqs. GJ-®. The result for the 
inverse structure factor matrix in reciprocal space can be 
written as 



(25) 



where first term comes from the ideal and correlation 
contributions to the free energy and the second term from 
the mean field electrostatics. The first term is in detail 



Tij 1 = Sij/Pi + p m Z 2 TT 2 l B a 3 hi((TK, erg) Ay 

- Z 2 irl%(xh2(<TK,aq)A'! j 
where the functions hi^(x — o-n,y — aq) are 

1_ (x 3 (x + 2f) ' 2 ~ {x{x + 2) 2 Y 
and the matrices are 



(26) 



A' = A' , 

mm m± 



A" — A" 



0. 

0, 



A'_. 



±± 



A" 



1, 
1. 



(27) 



(28) 



The y-dependence (y = aq) in Eq. I|27|l arises from the 
Fourier transform of the weight function of Eq. Note 
that the point model alluded to in section ITT1 corresponds 
to the limit a — > in Eqs. I|27l) . In this limit, the theory 
becomes ill-defined since Sij (q) does not have the correct 
limiting behaviour as q — > oo. This was the original 
technical reason for introducing the smoothing kernel. 

For any given state point and value of q, Eqs. I|25[l- 
PJl define S^ 1 which can be inverted numerically to 
find all components of the structure factor matrix. A 
partial solution can be obtained analytically in terms of 
the subsidiary matrix Ty, 



Si 



47r/ B J2ki z kZiT tk Tji 
q 2 + 4irl B J2m z kZiTki 



(29) 



From this one can readily prove that Sij exactly satisfies 
the Stillinger-Lovett moment conditions in Eqs. and 
J2U above. 

Another result follows from the dominance of the ideal 
contribution over the correlation contribution at low den- 
sities. In the limit pi — > one finds — ► pidij and 



Si- 



Pih 



47rl B ZiZjPiPj 
+ E fc z\pk 



(30) 



This is in fact exactly in accordance with the Debye- 
Hiickel limiting law at low densities. To see this, note 
that A = (47t^b Efc z fcPfe) -1 ^ 2 1S the Debye screening 
length defined to include all ionic species. Thus in 
real space, Eqs. ifTTI) and l(3TJI) indicate that h%j — 
—ZiZj(l-B/r)e~ r l x , in correspondence with the Debye- 
Hiickcl limiting law. 

It is clear that the moment conditions and the Debye- 
Hiickel limiting law behaviour follow from the construc- 
tion of the DFT to include a mean-field contribution sep- 
arately from the correlation term. This construction is 
in turn motivated by the expected behaviour of the di- 
rect correlation functions Cy (r) at r — > oo, as Evans and 
Sluckin have described 0|- The form of the correlation 
term is unimportant, so long as it is regular both at q — > 
and pi — > 0. 

For the remaining part, I now focus on the macroion 
structure factor S mm - Note that the theory includes the 
macroion-macroion electrostatic interaction explicitly in 
the mean field term, and an additional indirect interac- 
tion in the correlation term. The computation of S mm 
reveals the combined effect of these macroion interactions 
on the macroion correlations. 

Typically S mm has a 'hole' in reciprocal space for 
qa < 1. This corresponds to the macroion electrostatic 
repulsions. Within the correlation hole though, there is 
additional structure. This becomes particularly impor- 
tant in the vicinity of the phase separation region. Two 
kinds of behaviour are possible: at higher salt concentra- 
tions S mm rises to a maximum as q — ► 0, or at lower salt 
concentrations S mm acquires a peak at some q* > 0. In 
the phase diagram, the two alternatives are separated by 
a (macroion) 'Lifshitz line' {6f|, defined to be the locus 
of points for which dS mm /d(q 2 )\ q= o = 0. Fig. Ufa) shows 
the two behaviours for a pair of typical state points above 
and below the Lifshitz line, and Fig. Etb) shows the Lif- 
shitz line superimposed on the bulk phase behaviour. 

Also shown in Fig. [3(b) is the spinodal line computed 
from the bulk free energy in Eq. © of section IIIII One 
can check that S mm (q = 0) diverges on this spinodal line; 
in fact all the q = components of the structure factor 
matrix diverge because the determinant of S^ 1 vanishes. 
For salt concentrations above the Lifshitz line, this di- 
vergence at q = can be accommodated within the gen- 
eral behaviour of the structure factor. Of course, state 
points within the binodal are metastable so the diver- 
gence is strictly only visible as the upper critical point is 
approached. The fact that the structure factors diverge 
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on the spinodal line is no coincidence, since thermody- 
namic consistency by the compressibility route is assured 
for a DFT 0. 

What happens at salt concentrations below the Lif- 
shitz line? Here, the peak in S mm at q* > is found 
to diverge before the bulk spinodal line is reached. The 
shaded area in Fig. EJb) shows the region where this oc- 
curs. A divergence at a non-zero wavevector is indicative 
of microphase separation In this case one would 

expect a charge-density-wave (CDW) phase to appear 
[39. l68l | . The shaded region extends below the binodal 
for bulk phase separation, so the CDW phase should be 
observable in this part of the phase diagram. In fact the 
CDW phase will be found whenever the lower critical 
point lies below the Lifshitz line. The general idea that 
a critical point in a charged system can be replaced by a 
CDW phase was advanced by Nabutovskii, Nemov and 
Peisakhovich jH^El- 

The location of the Lifshitz line depends on the pa- 
rameter a which sets the range of the smoothing kernel 
w(r) in Eq. QJ. If a < 0.40 the Lifshitz line moves up- 
wards past the upper critical point, which would then be 
expected to be replaced by a CDW phase too. On the 
other hand if a ^ 3.6, the Lifshitz line moves downwards 
past the lower critical point. These critical values of a 
only depend on the coefficient of q 2 in the expansion of 
the Fourier transform of w(r) about q = 0. 

The Lifshitz line discussed here pertains to the 
macroion structure factor. Although slightly different 
Lifshitz lines are expected for each component of the 
structure factor matrix, the locus of state points where 
the peak diverges (either on the spinodal or on the bound- 
ary of the CDW phase) should be the same for all com- 
ponents. 

Whilst the Lifshitz line line marks an obvious change in 
the behaviour of S mm , the cross-over from monotonic to 
damped oscillatory asymptotic decay of the correlation 
functions hij(r) is determined by Kirkwood or Fisher- 
Widom lines in the phase diagram 

ELEL The 

difference between these is rather subtle and 
one might loosely cover both possibilities by the phrase 
'Kirkwood-Fisher-Widom' (KFW) line. The importance 
of the KFW line lies in the fact that it also governs the 
asymptotic decay of the interface density profiles, which 
behave in the same way as hij 71}. Thus the calcula- 
tions reported in section IIVI above, which assume that 
there is no oscillatory behaviour in the density profiles, 
requires as a necessary minimum that the coexisting bulk 
densities both lie above the KFW line. The location of 
the KFW line is governed by the poles of Sij(q) in the 
complex q plane, which are either purely imaginary or 
occur as complex conjugate pairs, and are the same for 
all components of f?y If the pole nearest the real 

q-axis is purely imaginary, then monotonic decay is ex- 
pected; conversely if a pair of complex conjugate poles is 
nearest the real g-axis, then damped oscillatory decay is 
expected j?2 . Determination of the KFW line is a hard 
numerical problem and has not been attempted for the 



present DFT. However the presence of a peak in S mm (q) 
on the real g-axis at q = 0, or at q* > 0, ought to be 
indicative of whether the pole nearest the real g-axis is, 
or is not, purely imaginary. Thus the Lifshitz line should 
serve as a guide to the location of the KFW line. In sec- 
tion \N] therefore, care was taken to make sure that the 
coexisting bulk densities lie well above the Lifshitz line. 



VI. DISCUSSION 

The paper presents a density functional theory (DFT) 
for a macroion suspension. The excess free energy cor- 
responds to the macroion self energy evaluated using 
Debye-Hiickel theory. These approximations render the- 
ory tractable without losing the basic phenomenology 
which resembles that of other studies. The advantage 
of a DFT is that one can compute the interface structure 
and surface tension between coexisting phases. The re- 
sults are in accord with expectations from previous work 
PJ. In particular, the electrical structure of the interface 
gives rise to a junction potential analogous to the Don- 
nan potential across a semi-permeable membrane. This 
arises from an electric dipole moment density (per unit 
area of interface), which appears because charge neutral- 
ity is locally violated in the vicinity of the interface. The 
surface tension is found to be of the order /cb T/a 2 . 

Structure factors can be computed from the DFT. 
These are found to obey the Stillinger-Lovett moment 
conditions, although this is not a stringent test of the 
theory. The structure factors reveal an interesting phe- 
nomenon, namely that oscillatory behaviour can appear 
in the (direct) correlation functions, particularly at low 
ionic strength. Indeed there may be regions of mi- 
crophase separation in the vicinity of the critical points, 
corresponding to the appearance of a charge-density- 
wave (CDW) phases. This phenomenon is peculiar to 
asymmetric charged systems [3(|, and is strictly absent 
in symmetric systems such as the restricted primitive 
model. In this respect, the possibility of CDW phases 
is correlated with the appearance of the junction poten- 
tial, which is also strictly absent in symmetric systems 
|35| . Given the approximate nature of the DFT, only cer- 
tain aspects of the present analysis might be expected to 
survive in a full treatment. One of these is an upturn in 
macroion structure factor at small q, even in the absence 
of a true miscibility gap. This would reflect an increased 
osmotic compressibility in this region of the phase dia- 
gram. Another expectation is the possible appearance of 
the CDW phases, although it might be difficult to disen- 
tangle these from the ordered (crystal) phases that are 
expected for a macroion suspension at sufficiently strong 
electrostatic coupling. 

The macroion self energy depends on the local ionic 
strength, but on both physical and technical grounds it 
is found necessary to introduce the notion of smooth- 
ing or smearing — the dependency should be on the ionic 
strength averaged over the vicinity of the macroion. Here 
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a completely phenomenological approach has been taken 
to construct the details of the DFT. Other choices could 
be made, or indeed more rigor could be introduced, such 
as additional requirements for internal consistency |74| . 
Tests indicate though that the general phenomenology 
(electrical structure of interface, gross behaviour of struc- 
ture factors) is found to be insensitive to the details of 
the model at this point. 
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APPENDIX A: CORRECTION TO REF. 

Chan has remarked that an excluded volume con- 
tribution was omitted in the theory of Ref. @- This 
appendix describes the missing term. The error occurs 
m going from Eq. (3) to Eq. (7) of Ref. @ where the 
omitted contribution arises from the fact that h m ±(r) = 
g m ±{r) — 1 = —1 for r < a/2. In terms of the microion- 
macroion interaction energy, i? ms / (V k^T) , the omitted 
contribution is 

f Zl 
Pm I d 3 r — - [p + h m+ (r) - p_h m _(r)\ 

J\r\<a/2 r 



cr/2 



-p m (p+ - P-) I 47rr dr ■ 

Jo r 



(Al) 



■nZ 2 l B p 2 m a 2 



(using p+ - p- = -Zp n 



This contribution is a positive, increasing function of 
p m , and has the tendency to stabilise the system against 
phase separation (because it is an athermal excluded vol- 
ume term, it passes unscathed through the thermody- 
namic integration step needed to calculate the contribu- 
tion to the free energy). If the calculations of Ref. [j| are 
repeated with this contribution included, it is found that 
the basic phenomenology is still the same, except that 
the miscibility gap in the (p m ,p s ) plane does not appear 
until somewhat larger values of ZI-q/cf. Fig. 8 shows the 
new results in comparison with those reported in Table 
II of Ref. |5|. The new calculation indicates that phase 
separation is observed in an even narrower window of pa- 
rameter space for which the Debye-Huckel linearisation 
approximation might be admissible than was found in 
the earlier work. This can be taken to indicate that the 
self-energy mechanism may not be sufficiently powerful 
to drive phase separation by itself, as discussed in the 
introduction. 
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FIG. 8: State diagram showing where a miscibility gap is 
found for the full theory of Ref. L 5| including the omitted 
term (solid line), compared to the original results (dashed 
line). The shaded region shows where Z > 4<t/?b, which 
is one possible criterion for the acceptability of the Debye- 
Hiickel approximation for the polarisation energy 



APPENDIX B: NUMERICAL APPROACH 

The task is to find density profiles pi(x) which min- 
imise the grand potential in Eq. (JTJ. The most accurate 
method is to solve the integral equations for the pro- 
files in Eq. Ijl3|l . However, this is hard. An alternative 
is to adopt a variational approach in which J7, or 7 in 
practice, is minimised with respect to parameters in trial 
functions which specify the density profiles |7^ . This is 
the approach that has been taken here. 

The ion density profiles have to satisfy a sum rule since 
the potential difference Aip = ijj(oo) — ip(— 00) is fixed by 
the coexisting bulk densities as described in section Hvl 
One can replace one of the ion density profiles by tp( x ) 
to ensure this sum rule is automatically satisfied. In the 
present case, a choice was made to use the set {p m , P+,4>} 
as a basis with p_ derived analytically from the Poisson 
equation, p_ = Zp m + p_ — (d 2 ip / dx 2 ) / (AitIb) ■ The first 
integral of the Poisson equation shows that one can ad- 
ditionally ensure global charge neutrality by making sure 
that dip/dx —* as |a;| — > 00. Once the pi are known, the 
average ionic strength 7J/ and the surface tension 7 are 
determined numerically by quadratures. 

To represent the basis set {p m ,p + ,^}, three copies of 
the function 



f(x;£±,{a}) = 



a+e 



-x/£_ 



a-a+ + a-e x ^+ + a + e- x /t- (Bl) 



are introduced. In this, the H r are Hermite functions, 
with £ = 2/(l/£_ + l/£+) used to scale the argument. 
Each copy of / is parametrised by the correlation lengths 
£± and amplitude set {a}, and has the properties that 
/ — > ±(1 — a±e Tx ^ ± ) as x — > ±00. One copy of / is 
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assigned to each member of {p m , p + , tp}, and is scaled 
and shifted to match the limiting values at |ar| — » oo, 
for example p m = p m (— oo)(l - /)/2 + p m (oo)(l + f)/2 
(for the electrostatic potential, one can set oo) = 
and i/j(oo) = Atp). The three copies of / have differ- 
ent amplitude sets {a} but share common values for £± 
since the asymptotic decay of the density profiles into 
the bulk phases is expected to be governed by a bulk 
correlation length — it is these values of £± that are re- 
ported in Fig.[|Jb). A finite set of N Hermite functions 
has been included in each copy of / to allow for an arbi- 



trary structure at the interface. In practice the minimi- 
sation problem is well behaved only if the density profiles 
smoothly interpolate between the bulk values, for which 
case typically N = 3-6 Hermite functions are needed to 
achieve convergence in 7 to an accuracy of the order 1%. 
At this point, the interface problem has been reduced 
to a multivariate minimisation over the three copies of 
the amplitude set {a} plus the correlation lengths £±. 
Numerical minimisation of 7 with respect to these pa- 
rameters is then undertaken by standard methods |76j . 
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